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ABSTRACT 

We present results of 2.5D numerical simulations of the emergence of 
sub-surface magnetic flux into the solar atmosphere, with emerging flux 
regions ranging from 10^® to 10^^ Mx, representing both ephemeral and 
active regions. We include the presence of neutral Hydrogen in the govern- 
ing equations, improve upon previous models by including the ionization 
in the equation of state, and use a more realistic convection zone model. 
We find that ionization and recombination of plasma during the rise of a 
convection zone flux tube reduces the rise speed of the tube's axis. The 
presence of neutral Hydrogen allows the effective flow of mass across field- 
lines, by the addition of a Pedersen resistivity to the generalized Ohm's 
law, which dissipates current perpendicular to the magnetic field. This 
causes an increase of up to 10% in the amount of magnetic in-plane flux 
supplied to the corona and a reduction of up to 89% in the amount of 
sub-surface plasma brought up into the corona. However, it also reduces 
the amount of free magnetic energy supplied to the corona, and thus does 
not positively affect the likelihood of creating unstable coronal structures. 

Subject headings: CMEs, Flux Emergence, MHD 



1. INTRODUCTION 



Coronal mass ejections (CMEs) and flares are the most energetic manifestations 
of solar activity. It is generally accepted that the energy required for these events 
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is stored as magnetic energy in the corona. This free magnetic energy is mos t hkely 
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tures have a significant component of the magnetic field parallel to the neutral line, 
compared to a potential configuration where the field is perpendicular to the neutral 
line. In the strong complex regions that are sources of fast C MEs, these filament 
channels form with the active region (jPeynman &: Martin 1995 ) and hence the shear 



field must emerge with the fiux. Given that shear formation is a fund amental driver of 



the e ruption in CME models such as t he magnetic breako ut model (lAntiochos et al. 



19991 ) and the fiux cancellation model (lAmari et al.ll2000l ). it is clear that fiux emer- 
gence must be explicitly included in such CME models for them to be physically 
rigorous. 



Most CME n iodels , such as those of lAntiochos et al.l (119991 ) . lAmari et al.l ( )2000h . 
Chen fc Shibatal J200oh . and IPan fc GibsoJ (120071 ) are constrained by the need to 
extend the simulation domain to at least a few solar radii. They therefore do not 
model the lower solar atmosphere. The lower boundary of these simulations has a 
typical density of 3 x 10~^^ g/cm^ and a typical plasma-/? < 0.2, where 



(3 



(1) 



P is the gas pressure, B is the magnetic field strength, and /xq is the permeability 
of free space. The ph otosphere, in contrast, has a density of around 3 x 10^'' g/cm^ 
( IVernazza et al.lll98ll ) and a /3 > 1 (IGary fc Alexanderlll999l : lGaryll200ll ). The supply 
of magnetic free energy to the corona in active regions is due to the emergence of sub- 
surface magnetic field structures, and so simulations of CME initiation must therefore 
address the critical question of whether and how newly emerging, sheared magnetic 
fiux can rise from its origins in the high /3 convection zone to the low corona where 
it is required to drive many CME models. 

The problem of fiux emergence has been studied extensively for a few decades. 
The difficulty in studying fiux emergence is that a simulation must couple the high /3, 
high density convection zone, where plasma pressure forces dominate, and the corona, 
which is low /3 and magnetically dominated. The basic idea of fiux emergence is that 
sub-surface fiux tubes are created by dynamo actions at the base of the convection 
zone. Indeed recent studi es suggest that deep convect ive layers control the pattern of 
large-scale solar activity (lArkhypov et al.ll201lU2012l ). These fiux tubes are assumed 
to have acquired sufficient twist to survive the convective motions as they rise to 
the surface. The subsequent expansion into the atmosphere and its interaction with 
both a field-free corona and pre-existing coronal fields has been stu died extens i vely i n 
recent years, and a comprehensive list can be found in the review of lArchontid ( 120081 ). 
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In lLeake et al.l (120101 ) we showed that emerging two and a half dimensional (2.5D) 
sub-surface flux tubes did not supply enough free energy to the corona (less than 20% 
of the magnetic energy was in the shear component of the field) and hence the re- 
sulting coronal magnetic field did not erupt. We concluded that in order to create an 
eruption a method is needed which allows the transfer of significant shear field into 
the corona by removing the dense material which is inhibiting the rise of the emerging 
structures. This can be achieved in two ways. One method is the draining of plasma 
along the axis of the emerging flux tube. Recent three dimensional (3D) simulations of 
flux emergence have shown that the drainage along strongly arched flux tube axes can 
aid th e rise of magnetic field into the corona (IHood et al.ll2009t iMacTaggart fc Hood 



2009bl ). Also, complex 3D motions associated with shear flows during flux emer- 



gence have helped incre a se the amount of free energy that emerges into the co rona 
^Manchester et al.l 120041 : krchontis fc Torokl boosi : IMacTaggart fc HoodI boOQah . A 
second method is the drainage of neutral plasma across field lines. Using a model 



for th e support of a prominence in a partially ionized solar atmosphere, [Gilbert et al 



( 120021 ) predicted the amount of vertical draining of neutral atoms from prominence 
structures, and also found evide nce of cross-field di ffusion of neutral material in solar 
filaments to support this model (IGilbert et al.ll2007l ). In this paper we will explore the 
drainage of neutral material across fieldlines during flux emergence, using a modified 
MHD model which includes the partial ionization of the solar atmosphere. 



In the simulations of iLeake et al.l (120101 ) and the majority of previous simula- 
tions, the numerical models use the magnetohydrodynamics (MHD) equations for 
a fully ionized plasma. For most of the convection zone and corona this is a valid 
approach. However, in the upper 3 Mm of the convection zone, as well as the pho- 
tosphere and the chromosphere, the temperature is low enough that the plasma is 
not fully ionized. It has been shown that the effects of ion-neutral co l lisions cre- 
ate a n anisotropic resistivity in the single-fluid equations ( jCowlingl 119571 : iBraginskii 



19651 ). This Pedersen resistivity acts only on currents perpendicular to the mag- 



netic field, and has been shown to be up to 12 orders of ma gnitude larger than 



the p arallel Spitzer, or Coulomb, resistivity in the chromosphere ( iKhodachenko et al. 



20041 ). The anisotropic dissipation of perpendicular currents by this Pedersen resistiv- 



i ty has been included in the study of the damping of MHD waves in the chromosphere 



(Ide Pontieu 



1999; 



Leake et al. 2005]). and flux emergence in both 2.5D and 3D simu- 



lations ( ILeake &: Arberll2006t lArber et al.ll2007l ). In those flux emergence simulations. 



it was shown that the Pedersen resistivity led to increased rates of flux emergence 
due to dissipation by ion-neutral collisions, and increased collisional heating. 



In the simulations presented here we im prove on two of t h e limi tatio ns of the par- 
tially ionized flux emergence simulations of iLeake fc Arberl (120061 ) and lArber et al. 
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(2007). The first is related to the equation of state used. In the simulations of 



Leake &: Arbeii ( 20061 ) and lArber et al.l ( 20071 ). the ionization fraction was calculated 
as a function of density and temperature, but the specific internal energy density 
(e) did not include a contributing term from the ionization fraction: The equation 
of state in those simulations did not include the change in internal energy density 
due to ionization and recombination. In the simulations presented here, this term is 
specifically included in the equation of state. As will be shown in the following sec- 
tion, this approach therefore also includes a correct description of a partially ionized 
convection zone which is adiabatically stratified. 

The second improvement is to use flux tubes which initially contain axial mag- 
netic fl ux comparable to obse rvati ons of real active r egions on the Sun. The simula- 
tions of lLeake fc Arberl (120061 ) and lArber et al.l (120071 ) use flux tubes in the convection 
zone with axial fluxes of less than 10^^ Mx. The axial flux in these tubes is two orders 
of magnitude less than that ty pically measured in a ct ive regio i is on the su rface of the 
Sun ( 10^^ Mx). More recently, Icheung et all toid ). iRempell toilh . andEng^taf 



(12012! ) have simulated the emergence of active region size flux tubes (> 10^1 Mx). 
However, those simulations do not include the effects of Pedersen resistivity on the 
flux emergence process. We model the emergence of large-scale tubes with axial 
fluxes of 1 0^^ Mx, as well a s the n iore typical small s cale flux tubes seen in t he sim- 
ulations of lLeake &: Arberl (120061 ). lArber et al.l (120071 ) . and lLeake et al.l (120101 ) . using 
our partially ionized plasma model. 



2. NUMERICAL METHOD 
2.1. Equations 

We solve the standard resistive MHD equations, modifled to include the effects 
of parti al ionization. These are solved numerically using the Lagrangian remap code 



Lare2d (lArber et al.ll200ll ). The equations are obtained by summing the equations 
for ions (i) electrons (e) and neutrals (n). Hence the total mass density (p), gas 
pressure (P) and internal speciflc energy density (e) are obtained by summing over 
the three species, e.g., p = YlaPa = 'Yl,a''^a.na where nia and Ua are the mass and 
number density, respectively, of each species (a = i,e,n). The average velocity, v, is 
deflned by pv = J^aPa^a,- The neutral fraction is deflned by 



= rin/iUn + Ui) . 



(2) 
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The equations are given below in Lagrangian form, using SI units: 

§ ^ -PV.V, (3) 

^ = _lvP+-jAB + g + lv.5, (4) 

Dt p p p 



B.V)v - B(V.v) - V A (r^j,,) - V A (r^pj^), and (5) 



Dt 

De P 

-j^ = -—^■^ + <iijSij + vj\\^ + Vpj±^, (6) 

where jy and jx are the parallel and perpendicular current vectors, respectively, and 
are defined as 

Jll = |g|2 ^ (") 

i BA(jAB) 
= ^"p • 

Here the current density is defined by j = V A B//io and /xq is the permeability of 
free space. The total gas density, pressure, and internal specific energy density are 
defined at the center of each numerical cell. The magnetic field, B, is defined at cell 
faces, and the velocity is defined at cell vertices. This staggered grid preserves V.B 
during the simulation. The gravitational acceleration is denoted by g and is set to the 
value of gravity at the mean solar surface (274 m/s^). S is the stress tensor which has 
components Sij = iy{Qj — |5jjV.v), with = + The viscosity u is set to 

3 X 10^ kg m~^s~^, and 6ij is the Kronecker delta function. The Coulomb resistivity, 
7], is given by 

, = !!^Mi±^^ (9) 

where mg and e are the mass and charge of the electron, respectively. The effective 
collisional frequencies for collisions of electrons with ions and neutrals, respectively, 
are u'^^ and z/g„. The Pedersen resistivity, rjp, is given by 

Vr = V + ^^- (10) 

an 

The quantity a„ is calculated using a„ = ■meneiy'^j^ + miUii'^^ where is the effec- 
tive collisional freq uency for collisions of ions with neutrals. We refer the reader to 



Leake et al.l (|2005[ ) for the method of calculating the effective collisional frequencies. 



For this partially ionized plasma, the total pressure, P, and the specific internal 
energy density, e, can be written as 

P = pkBT/fim, and (11) 
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TT + (l-a)-, 
1 rrii 



(12) 



respectively, where kB is Boltzmann's constant, 7 is 5/3, and Xj = 13.6 eV is the 
first ionizat ion energy of Hydrogen. Previous sira ulations of partially ionized flux 
emergence (iLeake fc Arberl l2006l : lArber et al.l 120071 ) did not include the second term 
in Equation f lT2|) . In this new model we allow for changes in energy to change the 
temperature and the ionization level simultaneously. The neutral fraction, is a 
function of temperature, T, itself and so Equation f|T2|) is solved implicitly for T at 
each time step. The reduced mass, fim, is 



/in 



mi/{2 - ^n) 



(13) 



Here = mfrrip, where rup is the mass of a proton, and = 1.25 is a pre-factor 
which is designed to include the effects of heavier elements and, as will be shown in the 
next section, will help reconcile our initial conditions with more realistic theoretical 
models of the solar convection zone. 

To calc ulate ^„ we use a simple model based on the modified Saha equation 
( lSahalll92ll ). which takes int o account th e fact that the chromosphere is not in local 
thermodynamic equilibrium ( lBrownlll973[). This equation can be s olved for the steady 
state solution of the ionization equation (lAthay &: Thomas! 1 1 96 ll ) to give 



where 

f{T) 
and 
h{T) 



m 

hiry 

{2TimekBT)~ 



T 

wTr 



exp 



■ exp( 



knT' 



(14) 
(15) 

(16) 



T 



AksT'TR 



Here Tr is the temperature of the photospheric radiation field, w is its dilution factor, 
and h is the Planck constant. Below the surface, Tji = T and w = 1 so that b{T) = 1. 
Above the surface, = 6420 K and w = 0.5. Given nf/nn-, the ratio of the number 
density of neutrals to number density of ions can be calculated from 



nn 




(17) 



and the neutral fraction, ^„ = nnjijin + rti), is 

r 



in 



1 + r 
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For the convection zone and lower solar atmosphere, the fully ionized value of 
the Spitzer/ Coulomb resistivity gives a magnetic Reynold's number = ULAh/c^t] 
(where U and L are typical measures of velocity and length) much larger than unity 
(> 10^). Thus typical numerical diffusion in MHD simulations of flux emergence 
excee ds the classical diffus ion, and hence = is often assumed. On the other 



hand iLeake fc Arberl (120061 ) showed that the value of the Pedersen resistivity gives 
an effective magnetic Reynold's number below unity in the chromosphere for flux 
emergence simulations, and that the increased diffusion due to the Pedersen resistivity 
is significantly larger than both the diffusion due to the Coulomb resistivity and 
typical numerical diffusion. 

The equations are solved in 2.5D: the simulation box is 2D, with x and y being 
independent variables and z being ignorable, but all three components of the vector 
variables are evolved. 



2.2. Background Stratification 

We require a background stratification which represents the solar convection zone 
and atmosphere. We first define a temperature profile which is chosen to resemble 
the real Sun. Then with a specific value of density (or pressure) at the surface, the 
hydrostatic equilibrium equation can be solved to give the density and pressure in 
the simulation domain. 

Above the surface (y > 0), a simple temperature profile is used 

T{y) = T,, + i^££^^M(tanh + 1), (19) 

/ Wtr 

where Tph = 5778 K, Tcor = 150 Tph, ycor = 3.75 Mm and Wtr = 0.75 Mm. This 
profile is desi gned to resembl e semi -empirical atmospheric models of the quiet Sun 



developed by IVernazza et al.l (Il98l[ ). Below the surface, an adiabatically stratified 



gas in hydrostatic equilibrium is used as a model for the atmosphere: The vertica l 



(y) temperature gradient is equal to the adiabatic temperature gradient (jStixll2004j ). 
In this way the atmosphere is marginally stable to convection. For a fully ionized 
plasma this yields a temperature profile in the convection zone of 



dyj a 7 



(20) 



Previous simulations of flux emergence which assume the plasma to be fully 
ionized typically use the neutral limit of the reduced mass of /x^ = rather than 
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/im = which is strictly the fully ionized limit. As we will show, this choice is 

made to ensure realistic coronal densities result from the solution to the hydrostatic 
equilibrium equation. 

Removing the assumption that the plasma is fully ionized yields a different result 
for the initial equilibrium. Rising and falling packets of plasma can now ionize and 
recombine, and do not expand and contract ideally. 

We want to derive an equation for , the adiabatic temperature gradient for 

a partially ionized plasma. This can be done by considering the adiabatic change of 
a volume V of partially ionized plasma. Using Equations ( IT2|) and ( ITSj) . and defining 
n = rii + Un, and ( = rii/n = 1 — ^n, we can express the energy in volume V as 

E^,.V^ ,„,nV ('-f^ + iV C-Sf^ + ex.) , (21) 



(7 - l)mi rrii J V (7 - 1) 

where N = nV is the total number of ions and neutrals and does not change. 

Adiabatic evolution assumes no heat exchange: dQ = 0. Using the first law of 
thermodynamics {dQ = dE + PdV), this implies 

dE = -PdV. (22) 

Differentiating Equation ( 12T1) gives 

dE = —^NkBdTil + C) + —^NkBTdC + NXidC- (23) 
7 — 1 7 — 1 

Writing the total pressure as P = ^^-bTXi+c) differentiate V to obtain 

dV f-dP dT dC \ 

' + — + T-h ) • (24) 

Inserting Equations (1251) and (1211) into Equation (1221) . and using P = ^^-bTXi+c) gjygg 



V \ P T 1+C 



where 



P "7-1 T + 1 + C' ^ ' 



^ + 7^- (26) 



7 - 1 ksT 

This is the condition for the adiabatic change for a partially ionized plasma, with 
arbitrary ionization state (. For practical purposes we would like to eliminate d( 
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from the equation, which we can do by rewriting the Saha equation (Eq. [T^ for 
y < in terms of (: 

= A- e-^^/'^-^, (27) 

\ — C p 

where A contains physical constants only. This can be differentiated to give 

(WT^^ - [-^ + ' ^^^^ 

or 

<K ^_f + »^^C(l^. (29) 



1 + C V ^ T J 
which can be inserted into Equation (1251) to give 



C(i-C) 



dT dP I 1 + 9 ^ 



T P \ ^ -i_ ^2 C(i-C) 

\ 7-1 ~^ 2 , 

Dividing by dy, and using ^ = —pg and P = pksT / gives 

dT\ fi^g f 1 + 



(30) 



(31) 



If C = or = 1, then we recover the usual definition for a fully neutral or fully 
ionized plasma, as expressed in Equation (!20|) . 

In summary, to construct a model atmosphere which is in hydrostatic equilibrium 
and whose temperature gradient below the photosphere {y = 0) matches Equation 
( 13T1) . we must perform the following procedure: 



1. For y > use the analytic profile T{y) from Equation f lT9|) . 

2. For the entire domain, initially set C = 1- Set T{y = 0) = 5778 K and p (y = 
0) = 3 .03 X 10~^ kg/m^, taken from the standard solar model presented in IStix 

3. For y < 0, numerically solve Equation f l3T]) with boundary condition T{y = 0), 
using the current C{y). This gives T{y < 0). 

4. For the whole domain, numerically solve the equation ^ = —pg with boundary 
condition p{y = 0), using the current T{y) and C{y), and the equation P = 
pkBT/pm- This gives p{y) for the whole domain. 
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5. Use the Saha equation (fT^ with the current T{y) and p{y) to get C{y)y and 
hence the reduced mass Pm{y) = i^^y) ^^e whole domain. 

6. Repeat steps 3-5 until the reduced mass Umiy) converges for all y. 

In order to show how this approach differs from the standard approach used in 
previous fully ionized simulations, we present the results of this process for three 
models. A summary of these simulations, along with simulations described in later 
sections, is shown in Table 1. Simulation uses the fully ionized plasma model with 
the ionized limit of /i^ = mi/ 2. Simulation 1 again uses the fully ionized plasma 
model but with the i ieutral lirait of / x^ = This is the case used by previous fully 
ionized simulations ( Archont"i^ boosi ) . We describe Simulation 2 later, and do not 



consider it at this point. Simulation 3 uses the partially ionized plasma model with 
an ionization which is allowed to change in space and time, based on the modified 
Saha equation. The results are shown in Figure [TJ Also shown are curves taken from a 
standard solar model which includes the transfer of heat by co nvection an d radiation 



in the partially ionized plasma of th e solar convection zon e (jStixl |2004J ). and a ID 



semi-empirical atmospheric model by lVernazza et al.l (Il98ll ). We can see from Figure 



[T] that the temperature of the partially ionized model (Simulation 3) best matches 
the standard solar model in the convection zone. The temperature of the two fully 
ionized models (solid and dotted lines) depart from the standard solar model curve 
in the convection zone, as the ionization cannot change in these two models. In all 
of the models we use frif = 1.25 in order to include the effect of heavier elements, 
and this ensures that at a height of y = the models have a reduced mass of 1.25mp 
which is the same as the standard solar model. However, only the partially ionized 
model is able to capture the variation in reduced mass with height in the convection 
zone and corona. 

While the temperature in the convection zone in the partially ionized model 
(Simulation 3) matches the standard solar model quite closely, there is some difference 
between the density and pressure profiles in the convection zone for the partially 
ionized model and the standard solar model. We expect that this difference is because 
the models here do not include the ionization of Helium, which is taken into account 
in the standard solar model. To test whether this difference was caused by the fact 
that the convection zone is not adiabatically stable throughout its depth, we also 
tested a model which had dT/dy = F{dT/dy)a where F = F{y) was a function 
w hich reprod uced the non-adiabatic temperature profile of the standard solar model 



of IStixl (j2004| ). This unstable convection zone profile did not produce any significant 



improvements in the density and pressure profiles. 
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Looking at the density and pressure at and above 2 Mm in Figure [H there is a 
large difference between the fully ionized model with the ionized limit of fim = '^j/2 
(Simulati on 0) and the other t wo models, Simulations 1 and 3. The atmospheric 
models of lVernazza et al.l (jl98ll ) give a density in the upper chromosphere/low corona 
(at a height of 2 Mm) of 10~^^ kg/m^, whereas Simulation has a density 4 orders 
of magnitude higher than this. Simulation was included here to show why previous 
authors who have performed fully ionized simulations use a reduced mass of /i^ = 
rather than /i^ = mi/2. Using the fully ionized limit of fim = mi/2 in a fully ionized 
simulation, which would seem the correct choice, actually creates large deviations 
from coronal densities. From now on Simulation is not used in this paper. 



EMERGENCE OF SMALL SCALE FLUX TUBES 



3.1. Initial conditions 



In order to cor npare our partially ionized simulations to previous small-scale flux 
emergence studies (iMagaral l200ll : iFanI l200ll : [Manchester et al.ll2004l ). we choose an 
initial flux tube with a very small axial flux compared to the flux of a real active 
region. In the next section we will address the emergence of flux tubes with active 
region size flux (10^^ Mx). We insert a cylindrical magnetic flux tube into the model 
convection zone at x = and y = yt = —1.8 Mm. The axial field, B^, and twist field, 
Bq, for the tube are given by 



B, = 5oe~"'/'^',and 
Bg = qrB^, 



(32) 
(33) 



respectively, where r = \/x^^^ljj'^^^^^ytj^ is the radial distance from the center of the 
tube, Bq is the axial field strength at the center (r = 0), a = 0.3 Mm is the radius of 
the tube, and the twist q = —1/a. For a field strength of Bq = 6000 G, this tube has 
an axial flux of 5 x 10^^ Mx. This tu be contains 25% less axial flux than our previous 
2.5D simulations (ILeake et al.ll2010l ). Our simulation domain is -3 to 42 Mm in y and 
-22.5 to 22.5 Mm in x. The numerical grid contains 640 x 640 cells, at a uniform 
resolution of 70 km. 

To investigate the effects of partial ionization, we perform three separate simula- 
tions. Simulation 1 is the fully ionized plasma model (with the neutral reduced mass 
limit of Jim. = mj), as describ ed in the previous section. To relate this work to the 
work of ILeake fc Arbed (120061 ) we also include the model used in that paper, and call 
it Simulation 2. In Simulation 2, the ionization is calculated as a function of p and 
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T using the Saha equation IHM . but the second term in Equation (fT2l) is ignored, so 
that ionization effects are not included in th e energy calculation. In addition, this 



means that the convection zone profile used in iLeake fc Arbeii (120061 ) does not include 
the variation of ionization in the equation for the adiabatic temperature gradient: it 
uses Equation fl20|) rather than Equation fl3Tl) . The initial background stratification 
for Simulation 2 is identical to Simulation 1, and so is not included in Figure [H Sim- 
ulation 3, as described in the previous section, uses the partially ionized model with 
the full equation of state, and includes ionization and recombination in the equation 
for the adiabatic temperature gradient. As shown in Figure [H Simulation 3 has a 
higher gas pressure at -1.8 Mm, where the flux tube is located, compared to Simula- 
tion 1. So for the same field strength, the flux tube will have a higher plasma /3 in 
the partially ionized simulation (3). We therefore vary the magnetic field so that the 
plasma /3 has the same value at the center of the t ube in all three simul ations. The 



rise speed of a buoyant flux tube scales as l/v^ (ILongcope et al.lll996l ). and so by 
matching /3 across simulations we ensure that the rise speeds are initially the same. 
For Simulations 1 and 2 we use Bo = 6000 G. For Simulation 3 we use Bq = 10900 G. 
The plasma (3, calculated as /3 = fioPoiVtube)/ B^{x,ytube), for the initial conditions of 
these three models is shown in Figure [2], Panel (a). 

In all these simulations we add a perturbation Pi{r) to the background plasma 
pressure Po{y) such that (VPi)^, = ( j A B)^, so that the tube is in radial force balance. 
Assuming that the flux tube is in thermal equilibrium with its surroundings makes 
the tube less dense than the surrounding plasma and initiates its buoyant rise to the 
surface. With this "isothermal" assumption, the density perturbation, pi(r), which 
is added to the background density po{y), is given by 

pi/po = Pi/Po ~ -1/2/3. (34) 
This perturbation for the three simulations is shown in Figure [2], Panel (b). 



3.2. Emergence from the convection zone into the atmosphere 



In the three simulations the flux tubes rise to the surface due to their initial 
buoyancy and undergo a degree of horizontal expansion as they meet the convectively 
stable photosphere. A contact layer is created when the tube's field meets the pho- 
tosphere. As more field builds up at the surface, t he layer becomes unstable to the 
magnetic Rayleigh Taylor instability ( lAchesonlll979l ). The criterion for this instability 
can be written as 

"-^^ ^;v> (35) 



- Hp—ln{B) > Hie 



1 + 
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(jArchontis et al.ll2004j) wher e Hp is the local pressure scale height and 5 is the super- 



adiabatic excess ( IStixl |2004J ). which is the double logarithmic temperature gradient 
{dlnT/dlnP) minus its adiabatic value. In these simulations, 5 is zero in the model 
convection zone, but negative in the model photosphere. The wavenumber k is the 
wavenumber parallel to the magnetic field in the horizontal plane, / is the wavenumber 
perpendicular to the field in the horizontal plane, and n is the wavenumber in the 
vertical direction. For the contact layer that is created, the left hand side in Equation 
(l35l) is positive, and acts to destabilize the layer. The last term on the right hand 



side is also positive and acts to stabilize the layer. 

Figure [3] shows the rise of the fiux tubes in the convection zone in the three 
simulations. The contour lines show fieldlines in the plane, given by constant intervals 
of where is defined as the vector potential for the 2D vector {B^-, By) = VAA^e^ 
with boundary condition ^^(oo) = 0. The outer fieldlines shown here are the contours 
of Az = OAmin{Az) where Az is always negative. The background color contour 
shows log(p X m^/kg). The tubes have risen to the surface with comparable speeds, 
which is to be expected as they started with the same plasma (3. The outer fieldlines 
are a different shape in Simulation 3, the partially ionized model simulation, than in 
Simulations 1 and 2. This is due to an increased horizontal expansion at the top of 
the tube in Simulation 3 relative to Simulations 1 and 2. This vertical gradient in 
horizontal expansion is caused by changes in ionization, which are not taken account 
in Simulations 1 and 2. 

Figure H] shows the height of the center and edge of the tube in all three simu- 
lations. The center of the fiux tube is defined as the location of the local extrema 
of Az, as Az is initially a negative monotonically increasing function of radius for 
our twisted fiux tubes. The edge of the fiux tube is defined as the intersection of 
Az = 0.1min(y42) with the a; = line. In Simulation 3 (dot-dashed line), the fiux 
tube's center stops 0.5 Mm below the surface. On the other hand, in the fully ionized 
simulation (1) and the partially ionized simulation (2), the tube's center rises toy = 
(the surface). Figure HI Panel (b) shows the rise speed in the convection zone in all 
three simulations. All three simulations show an identical rise speed initially, as they 
have the same plasma /3. The rise speed in Simulation 3 slows down relative to that 
of Simulations 1 and 2 due to changes in ionization during the rise. 

At around 1000-1020 s, all three simulations show the onset of the magnetic 
buoyancy instability. Figure O shows the stabilizing {—^(36) and the destabilizing 
{—Hp-^ln{B)) terms in Equation (135!) for the three simulation at two different times. 
As magnetic field reaches the photosphere, the local plasma-/? decreases, which re- 
duces the stabilizing term —^f36, until it becomes less than the destabilizing term 
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—Hp-^ln{B). At this point, the instabihty starts to develop, and magnetic field 
emerges into the atmosphere. Simulation 2 shows a much smaller gradient in mag- 
netic field than the other two simulations, i ndicated by a lower m agnitude of the 



destabilizing term in Figure El Panel (e). In iLeake fc Arberl ([20061) we showed that 
the magnetic field was able to slip through the plasma due to the Pedersen resistivity, 
which dissipated currents perpendicular to the field. In Simulation 2 we see the effect 
of this as a reduction in the vertical gradient of magnetic field in a layer above a height 
of 0.6 Mm above the surface. We do not see this in Simulation 1 as the model does 
not include Pedersen resistivity in that simulation. At this point in time, we do not 
see this in Simulation 3, as the rise of the tube in the convection zone in Simulation 
3 was slowed due to the ionization/recombination effects, but this reduction in the 
gradient in magnetic field is seen later in the simulation. 

Figure [6] shows the expansion of the magnetic field into the corona due to the 
magnetized Rayleigh Taylor instability, and Figure [7] shows the magnetic Reynolds 
number 

Rm ~ ^^^A^ (36) 

for the two partially ionized simulations (at time 1359 s for Simulation 2, and at 
time 1724 s for Simulation 3). The fully ionized simulation has a theoretically infinite 
Reynold's number (due to rip = 0). At around 0.3 to 0.8 Mm above the surface, 
the two partially ionized simulations show a value of Rm < 1- This means that the 
field's motion is dominated by diffusion rather than advection, and the dissipation 
increases the emergence of the fi eld into the atmosphe re. This result is in agreement 



with the previous simulations in lLeake fc Arberl (120061 ). even though in Simulation 3 



we have a different equation for the specific energy density e which allows for changes 
in ionization. 

As shown in Figure IH the upper edge of the tube expands further and faster 
into the corona in Simulations 2 and 3, as the increased dissipation due to ion-neutral 
collisions allows the field to diffuse through the plasma without having to lift material 
up. In the fully ionized simulation, the field must lift more plasma up with it, which 
slows its rise. Figure [H] shows the amount of in-plane fiux that has emerged into the 
corona for each simulation, normalized to the initial amount of in-plane fiux in each 
flux tube. This in-plane flux is calculated by 

<l>(t) = [max{A^{x,y)) - min{A^{x,y))]y>y^ (37) 

where yi = 1.2 Mm. The flux $(t) is normalized to \min{Az)\ at t=0, so that it 
represents the fraction of initial in-plane flux that emerges above the height 1.2 Mm. 
All three simulations emerge above 60% of the initial in-plane flux above this height. 
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Due to the Pedersen dissipation, the two partially ionized simulations have peaks in 
the in-plane flux above 1.2 Mm higher than Simulation 1, the fully ionized simulation. 

In the fully ionized model, the expansion of the plasma due to the emerging 
field is associated with a significant amount of cooling, and the plasma temperature 
drops to around 100 K within the emerging flux region, which can be seen in Figure 
M The transition region is pushed outward from 2 Mm to 14 Mm. The ion-neutral 
collision effects in the two partially ionized models (Simulations 2 and 3) are able to 
counteract this expansive cooling, and so this amount of cooling is not seen in the 
two partially ionized simulations, although some cooling to 1000 K is still seen and 
the transition region is still pushed outwards. There are two effects which contribute 
to this result. The first is collisional ion- neutral heating. The second effect is the 
slippage of mass through the field, due to the Pedersen resistivity, which allows the 
field to expand without expanding and cooling the plasma. These two effects ensure 
photospheric/chromospheric tempe ratures stay above 100 K in the partially ionized 



simulations, as originally shown in iLeake fc Arberl ( 120061 ). This result gives further 



evidence that ion-neutral collisions are an important source of heating in emerging 



active regions. Radiative-MHD numerical investigations by iLeenaarts et al.l (1201 if ) 
show that, in 2D simulations, unrealistic cooling occurs which can only be counter- 
acted by Joule heating. We have shown here that the additional Joule heating due 
to ion- neutral collisions is one way to achieve this necessary effect. 

We now investigate the effect of the increased dissipation due to ion-neutral 
collisions on the amount of mass and free energy supplied to the corona during flux 
emergence. In Simulation 3, the density is higher at the initial tube location, as 
shown in Figure (H and hence the flux tube contains more plasma than in the other 
two simulations (1 and 2). We therefore calculate the change in mass in the corona, 
above a certain height during the flux emergence, normalized to the amount of mass 
initially in the flux tube. This diagnostic effectively estimates the fraction of the flux 
tube mass which is supplied to the corona. Figure [10] shows this diagnostic above 
1.2 Mm in the three simulations. This shows that at t=2900 s, the fully ionized 
simulation (1) initially has lifted up approximately 7.7 times more of the tube's mass 
above 1.2 Mm than Simulation 3, but only 28% more than Simulation 2. 

Figure [TT] shows the quantity jjlj^l, which is equivalent to sin(6'), where 6 is the 
angle between the current density and magnetic fleld, for the three simulations, during 
the post-emergence stage. A value of means that all the currents are aligned with 
the fleld, i.e., the fleld is force-free. This shows that the angle between the fleld and 
the current density is lower inside the expanding tube in Simulation 3 compared with 
Simulation 1. The Pedersen dissipation acts on the cross-fleld currents, thus reducing 
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■|4j|^. This reduction in the Lorentz forcing component of the current density in part 
explains why so much less mass is supplied to the corona in Simulation 3, relative 
to Simulation 1. Simulation 2 also shows lower values of j4|j^ inside the emerging 
active region, relative to Simulation 1, again, explaining in part the reduction of mass 
supplied to the corona. 

It is worth noting here that the Pedersen resistivity in the partially ionized model 
only effectively dissipates perpendicular currents inside the emerging flux structure. 
The edge of the flux tube is a current sheet, with a large perpendicular current. How- 
ever, for this current sheet, where the magnetic fleld decreases to zero, the Pedersen 
resistivity falls to zero as the current increas es, and so there is limited dissipation of 



the current sheet. In fact lArber et al.l (120091 ) showed that for current sheets with no 



guide fleld, the Pedersen resistivity acts to sharpen currents, rather than dissipate 
them. 

There are two effects of the ion-neutral collisions working in parallel in the two 
partially ionized models during the emergence. Firstly, due to the Pedersen resistivity, 
the mass can 'slip' through the emerging magnetic fleld, and so less mass is lifted per 
unit amount of concave up flux in the partially ionized model relative to the fully 
ionized model. Secondly, cross-fleld currents are reduced by the Pedersen resistivity, 
and so the magnetic fleld cannot lift up as much mass in the partially ionized models 
as in the fully ionized model. For the parameters chosen for these simulation, we flnd 
that the result of these two effects is the reduction of the amount of mass supplied 
to the corona in both partially ionized models (2 and 3), relative to the fully ionized 
simulation (1). 

Figure [10] also shows the amount of energy in the shear component of the fleld 
{Bz) normalized to the total amount of energy in the fleld. 



p=22.5Mm ry=42Mm tj2j^j„ 
rp _ Jx=-22.5Mm Jyi -P^MXat/ 

^shear — r.x=22.5Mni ry=42Mm |t3|2 j j ^ 
Jx=-22.5Mm Jyi \ri\ dXdy 



where yi = 1.2 Mm. In the fully ionized simulation (Simulation 1 - solid line), the 
shear fleld {B^) is coupled to the mass, and we see more shear fleld supplied to the 
corona than in the other two simulations. At t=2900 s. Simulation 1 has a value 
of Eshear 10% higher than Simulation 2, and 25% higher than Simulation 3. We 
postulate that Simulation 1 can emerge more shear flux, even though this means 
emerging more mass, because it has more Lorentz force. 



Comparing Simulations 2 and 3 shows that including ionization effects in the 
equation of state signiflcantly affects the emergence of flux. There is more horizontal 
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expansion during the emergence process, both in the convection zone and the corona 
when the full equation of state is used. By 2900 s, Simulation 2 has emerged 6 times 
the normalized mass into the corona that Simulation 3 has. Also, Simulation 2 has 
emerged 15% more shear (in terms of E^^he ar) into the corona th an Simulation 3. 



We also note that the conclusions made in iLeake &: Arberl ( 120061 ). which used the 
model of Simulation 2, are repeated here with Simulation 3, which uses a more self- 
consistent model for the partially ionized plasma, in that it includes instantaneous 
ionization/recombination in the equation of state. These results include the increase 
in the amount of in-plane flux injected into the corona (Simulations 2 and 3 emerge 
more in- plane flux than Simulation 1), and the additional Joule heating which is vital 
for maintaining a chromospheric temperature above 1000 K. 

In summary, we have performed three simulations of an emerging flux tube with 
a single value of /3. Simulation 1 was a fully ionized model. Simulation 2 was a 
partially ionized model without the full equation of state, and Simulation 3 was a 
partially ionized model with the full equation of state, and a correct treatment of the 
adiabatic convection zone. We have shown that if a fully ionized model is used, the 
neutral limit for firn = should be used to get the correct density in the corona. 

These results show that although the presence of neutral Hydrogen can reduce 
the amount of mass supplied to the corona, which in principle can increase the likeli- 
hood of eruption, the amount of shear supplied to the corona in the partially ionized 
simulations is actually less than in the fully ionized simulations. Hence the likelihood 
of eruption is less in the partially ionized simulations. In all three simulations, the 
amount of magnetic energy in the shear component of the field is less than 20% of 
the total magnetic energ y in the field , and t hus the coronal field is non-eruptive, just 



as in the simulations of iLeake et al.l (|2010[ ). For these flux tubes, the 'slippage' of 



field through the nearly neutral lower atmosphere does not positively increase the 
likelihood of eruption. 



4. EMERGENCE OF LARGE SCALE FLUX TUBES 

The results of the previous section showed that including the ionization in the 
equation of state of the partially ionized model gave quantitatively different results 
from the model without the ionization terms. We performed t he comparison in orde r 



to put the results in context with the previous simulations of iLeake fc Arberl ( 120061 ). 
We know that the model used in Simulation 3, which includes the ionization in the 
equation of state, is the more self-consistent model of the two partially ionized models, 
and so we drop the model used in Simulation 2, which does not include the ionization 
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term in the equation of state. We study now the fully ionized model and the partially 
ionized model in terms of the emergence of active region size magnetic flux tubes. 

The axial flux in the tube in the previous section's simulations was less than 10^^ 
A/r,, „rV,;„v, ^,,„v, r,^^^^^^ than the 10^^ Mx size of a typical active region sunspot 



(Priest 


2003; 


Zirin 


1998) 



values to better test the effect of partial ionization on emerging active regions. 

To create an active- region scale flux tube, we increase the tube radius to 1.5 Mm 
and place the tube at a depth of -4.95 Mm. The domain is increased to cover a range 
-56.25 to 56.25 Mm in x and -15 to 97.5 Mm in y. The numerical resolution is kept the 
same as before at 70 km. The results of the previous section showed that flux tubes 
with initially comparable /9 in the flux tube start to rise with comparable speeds. In 
this section we run four simulations, which are grouped into two pairs of comparable- 
/3 simulations. The first pair consists of a fully ionized model (Simulation 4), and a 
partially ionized model (Simulation 5). Both simulations start with a fiux tube which 
has /3 = 4. The second pair consists again of one fully ionized model (Simulation 6) 
and one partially ionized model (Simulation 7), both having a fiux tube /3 of 40. The 
simulations are briefiy described in Table 1. The axial fiux in the tube in Simulation 
5 is 1.4 X 10^^ Mx, and so is comparable with observed active region sunspots. The 
two different values of /3 of 4 and 40 are chosen to cover a range of active region 
formation time. The simulations from the previous section show that fiux tubes with 
an initial plasma f3 of 2 emerge in less than 30 minutes, and so increasing the plasma 
/3 in the tube to 40, we expect emergence on a timescale of approximately 2 hours, 
which is somewhat closer to the general formation time of active regions on the Sun 
(12-48 hours). Figure [12] shows the initial (3 profiles and density perturbations, pi/po, 
for these four simulations. 

Figure [13] shows the in-plane magnetic field as regular contours of Az, and log of 
density as a color contour, at two different times in Simulations 4 and 5 (the /3 = 4 
simulations). The first time (left panels) is when the tube has reached the surface, 
and built up enough magnetic field to trigger the magnetic buoyancy instability. The 
second time (right panels) is at a later stage, when the outer 10% of the in-plane 
fiux has reached approximately 10 Mm above the surface. Figure [H] shows the same 
phases of the emergence process but for the two simulations 6 and 7, where /3 = 40. 
The high /3 in these simulations means that the initial rise-speed is slow compared 
to the low P simulations 4 and 5, and it takes longer to build up enough field at the 
surface to trigger the instability which drives field into the corona. 

Figure [TS] Panel (a) shows the height of the tube center and tube edge for the 
large scale simulations. As in the small scale simulations, the tube centers are con- 
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fined to just below the surface, wliile tlie edges expand in to tlie field-free corona. 
Figure [151 Panel (b) shows the rise speeds in the convection zone. Firstly, we note 
that these simulations confirm that the convection zone rise speed scales with 1 / a/^ 



( iLongcope et al.lll996l ): Taking the two fully ionized simulations (4 and 6), and cal- 
culating the ratio of the maximum rise speeds in the convection zone gives a value 
of 3.3, which is approximately equal to the inverse of the ratio of a/5 in those two 
simulations (a/40/4). The same applies to the rise speeds in the two partially ionized 
simulations (5 and 7). 

Comparing Simulation 4 to Simulation 5 (the two (3 = 4 models, solid and dot- 
dash line respectively), the rise speed peaks at a value of 2.3 km/s in the partially 
ionized model (Simulation 5), and a value of 2.9 km/s in the fully ionized model (Sim- 
ulation 4). This difference in rise speeds also occurred in the small scale simulations, 
and is a consequence of the ionization and recombination during the rise of the tube 
in the partially ionized convection zone. This result is repeated in the /3 = 40 simu- 
lations (Simulation 6: double-dot-dash line, and Simulation 7: long-dash line), with 
the fully ionized simulation obtaining a higher rise speed than the partially ionized 
simulation. This may explain why the flux tube axis settles lower in the convection 
zone in Simulation 7 compared to Simulation 6. This effect may be important for the 
formation of coronal structures in 3D, as it is thought that the evolution of the flux 
tube axis is important for the formation of sheared structures in 3D. 

The amount of normalized in-plane flux ($) supplied to the corona above 1.2 Mm 
is shown in Figure [161 Comparing $(t) for the low and high /3 simulations, we see 
that more of the original in-plane flux emerges in the low P simulations (4 and 5) than 
in the high /3 simulations (6 and 7). Still there is almost 50% of the original in-plane 
flux supplied to the corona in the high /3 case, which is comparable to the amount 
of flux which emerges above 1.2 Mm in the small-scale simulations of the previous 
section. Comparing the fully ionized simulations to the partially ionized simulations 
(4 vs. 5 and 6 vs. 7 ) it is clear that the flux emerges earlier in the partially ionized 
simulations (5 and 7), compared to the fully ionized simulations (4 and 6). 

The resulting temperature profiles are shown in Figure [T71 As in the small scale 
simulations, the partially ionized simulations do not see the drastic cooling below 
1000 K in the chromosphere that is seen in the fully ionized simulations. This is 
a result of both the collisional heating, and the dissipation by Pedersen resistivity 
reducing the expansion of the plasma as the field expands, a result which we found in 
the smaller-scale simulations of the previous section. We have therefore shown that 
this result applies to the larger active region and slower emergence of these large scale 
simulations (4 through 7). 



- 20 - 



Figure [18] shows the normahzed change in mass above 1.2 Mm in the four large 
scale simulations. Comparing the two /3 = 4 simulations (4 and 5), we find that the 
fully ionized simulation supplies approximately 10 times more of the fiux tube mass 
to the corona than the partially ionized simulation, a result which is repeated for the 
two /3 = 40 simulations (6 and 7). This, again is a direct consequence of the Pedersen 
dissipation in the partially ionized models. 

It is worth comparing the amount of fiux tube mass supplied to the corona in 
these large scale simulations to that in the small scale simulations. The two large 
scale (3 = 4 simulations supply approximately the same amount of fiux tube mass 
above 1.2 Mm (10~^ for the fully ionized model and 10~^ for the partially ionized) 
as the small scale simulations, which have /3 ~ 2. The /3 = 40 large scale simulations 
supply more than an order of magnitude less fiux tube mass to the corona than the 
/9 = 4 simulations (5 x 10^^ for the fully ionized and 4 x 10^*^ for the partially ionized 
model). So in general the amount of fiux tube mass supplied above 1.2 Mm during 
flux emergence is inversely proportional to the initial /3 in the fiux tubes at t=0, and 
hence is dependent on the strength of the magnetic field in the tubes, as one would 
expect. 

Also shown in Figure [18] is E shear above 1.2 Mm. For the /3 = 4 simulations, 
the fully ionized simulation supplies 15% more free energy at t= 2900 s than the 
partially ionized simulation. There is no significant difference between the two /3 = 40 
simulations. In these large scale simulations, less than 25% of the magnetic energy 
in the corona is in the shea r component , whic h broadly agrees with the results of the 



fully ionized simulations of iLeake et al.l ( 120101 ). Although the free energy does exceed 



20% of the magnetic energy in the /9 = 4 simulations, it quickly falls below this value. 
The /3 = 40 simulations stay below 20% for the duration of the simulations. The 
resulting coronal configurations in these 2.5D simulations are non-eruptive, despite 
the decrease in mass supplied to the corona due to the presence of neutral material 
in the lower solar atmosphere. 

The successful emergence of fiux tubes with active region values of magnetic fiux 
(10^^ Mx) has been achieved, and been shown to be qualitatively similar to the emer- 
gence of smaller scale fiux tubes, which have been studied extensively. Furthermore, 
the rise speed of these large scale fiux tubes is affected by the partially ionized region 
of the solar atmosphere, as is the amount of mass and shear supplied to the corona, 
even though the initial tube radius (1.5 Mm) is comparable to the size of the region 
where the Pedersen resistivity is important (from the surface to about 1 Mm above 
the surface). 
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5. CONCLUSIONS 



We have investigated the effect of neutral Hydrogen in the lower solar atmosphere 
on the emergence of sub-surface magnetic flux into the solar atmosphere, using 2.5D 
MHD models modified to include a variable ionization state, and the anisotropy 
in Ohm's law caused by ion- neutral co llisions. Previous simulations using a fully 
ionized MHD model (ILeake et al.ll2010l ) have shown that the amount of free energy 
is insufficient to drive a CME when flux tubes emerge in this 2.5D cartesian setup. 
The presence of too much mass and the lack of a strong enough shear field yielded 
no further expansion of the flux rope structure, which was ultimately constrained by 
overlying field. Typically the maximum amount of energy in the shear field was 10- 
15% of the magnetic energy. The main aim of this paper was to investigate the effect 
of the partially ionized atmosphere on the likelihood of creating eruptive magnetic 
field from emerging convection zone flux tube. 

The single-fluid MHD equations were modified to include the presence of neu- 
tral Hydrogen. This led to a modified induction equation, where perpendicular cur- 
rents were dissipated by ion-neutral collisions and parallel currents were dissipated 
by electron-ion and electron-neutral collisions. It also led to a source term related 
to ionization state in the equation f or the specific internal energy density. Pre- 



vious simulations (ILeake et al.l l2005l : iLeake fc Arberl l2006l : lArber et al.l 120071 ) with 



these effects used a simple equation of state which did not take into account ioniza- 
tion/recombination, and thus used a model for an adiabatic convection zone which 
did not include partial ionization. In this paper we included the change in ionization 
in our equation of state, and therefore a more realistic convection zone profile. 

In the partially ionized simulations, the dissipation due to the Pedersen resistivity 
and the associated coUisional heating was concentrated in the lower atmosphere, below 
the transition region. In the fully ionized simulations, the rapid expansion of the 
emerging field was associated with a rapid expansion and cooling of the plasma in the 
emerging region, which resulted in temperatures below 1000 K, which we consider to 
be too low for the real Sun. In the partially ionized simulations, we did not see such 
cooling below 1000 K. There are two effects which contribute to this. The first is the 
fact that as the field emerges into the corona in the partially ionized simulations, the 
Pedersen resistivity allows the field to 'slip' through the plasma, and this reduces the 
amount of expansion and cooling seen in the emerging region. The second is that 
any cooling will be countered by the collision dissipation. This result gives further 
evidence that ion-neutral collisions are an important source of heat ing in emerging 



active regions, as was suggest in the radiative-MHD simulations of iLeenaarts et al 
J201lh . 
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The result of the increased dissipation in the lower atmosphere in the partially 
ionized simulations led to a number of effects. Firstly, for the small scale simulations 
and the /3 = 40 large scale simulations, the amount of in-plane flux supplied to the 
corona increased compared to the fully ionized models. The Pedersen dissipation in 
the partially ionized simulations led to a decrease in the Lorentz forcing component 
of the current density inside emerging active regions, and a reduction in the mass 
lifted into the corona by the emerging flux. In the fully ionized simulations, there was 
more of the flux tube mass supplied to the corona. Coupled to this was an increase 
in shear field supplied to the corona in the fully ionized simulation. As a result, the 
likehhood of eruption was not increased by the effects of neutral Hydrogen in the 
lower atmosphere, but instead reduced. 

We performed studies with both small scale (radius of 0.3 Mm and fluxes of 
5 - 10 X 10^^ Mx), and large scale (radius of 1.5 Mm and fluxes of 10^^ - 10^^ Mx) 
flux tubes. The larger flux tubes resulted in active regions of size 20 Mm, based on 
the separation of the two intersections of the 10% flux contour with the surface y = 0. 
For the smaller flux tubes the active region size was typically 10 Mm. Interestingly, 
increasing the initial flux tube radius by a factor of 5 only resulted in an increase in 
the resultant active region size by a factor of 2. 

Results of our small scale fiux emergence simulations showed that by including 
the ionization in the equation of state, the rise speed was cfTcctivcly reduced by 
ionization and recombination in the upper convection zone, leaving the axis of the tube 
lower down in the convection zone after the emergence. Also in the simulations where 
the ionization was included in the equation of state, the amount of mass and shear 
flux supplied to the corona was decreased relative to the partially ionized simulation 
which did not include ionization in the equation of state. 

For the larger scale simulations we used two different vales for (3 in the flux tube. 
The timescale of emergence depended on the (3 value. For ^ — 4 the flux tube edges 
took 50 minutes to reach a height of 35 Mm, whereas for P — AO the emergence took 
100 minutes to reach the same height. This delay was almost entirely accounted for by 
the slower rise speed in the convection zone, which scaled with a/1//3. In all the low 
and high f3 simulations, more than 55% of the initial in-plane flux emerges into the 
corona. Even though the flux tube's radius in the large scale simulation was 10 times 
larger than the photospheric scale height of 150 km, the emergence was qualitatively 
similar to the emergence of smaller flux tubes. 

For these 2.5D simulations, the presence of neutral Hydrogen in the lower atmo- 
sphere, and the associated 'slippage' of emerging magnetic field (or equivalently the 
associated motion of plasma across fieldlines) did not increase the amount of shear 
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flux supplied to the corona, and therefore we conclude that the effects of neutral Hy- 
drogen do not increase the likelihood of eruption. We also conclude that 3D plasma 
motions along the axis of the flux tubes are more likely to destabilize coronal mag- 
netic field in emerging active regions. However, the stability of fiux ropes formed 
during fiux emergence, and the subsequent likelihood of eruption, will be dependent 
on the currents associated with the emerging fields, and therefore should be strongly 
dependent on the nature of the current dissipation mechanism. We propose in future 
work to investigate both 3D effects and the effects of partial ionization. Note that 
these conclusions address how the presence of neutral Hydrogen affects the supply of 
free energy, or sheared field, into the corona by the emergence of magnetic fiux into a 
field-free corona. Hence the only possible source of free energy in these simulations is 
the newly emerging flux. In the case where a pre-existing sheared structure is already 
formed in the corona, the emergence of new fiux can play the role of destabilizing the 
magnetic field and causing an eruption. 
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Name 


Model 


Reduced 
mass /im 


Eq. of state 


CZ profile 

(f),. 


Tube radius 


BoikG) 


Ptube 


Flux 

xlO^* Mx 





FIP 


mfmp/2 




Eq.(l20]) 










1 


PIP 


mfirtp 




Eq.(l20]) 


0.3 Mm 


6 


2 


5.4 


2 


PIP 


Eq.([l3D 


No 2nd term 


Eq.([20]) 


0.3 Mm 


6 


2 


5.4 


3 


PIP 


Eq. W) 


2nd term 


Eq.([3l]) 


0.3 Mm 


10.9 


2 


9.8 


4 


FIP 


mfiTLp 




Eq.([20D 


1.5 Mm 


15 


4 


340. 


5 


PIP 


Eq.(Il3D 


2nd term 


Eq.([3l]) 


1.5 Mm 


61.2 


4 


1380. 


6 


FIP 


rrifmp 




Eq.(l20]) 


1.5 Mm 


0.49 


40 


11. 


7 


PIP 


Eq.(Il3]) 


2nd term 


Eq.(l3l]) 


1.5 Mm 


19.2 


40 


430. 



Table 1: Simulations used in this paper. Simulation is used only to highlight the 
choice of the reduced mass /im and is not run beyond t = s. Sim ulation 2 has 



the same background stratification as the model used in lLeake fc Arberl (120061 ). FIP 
stands for 'fully ionized plasma', and PIP for 'partially ionized plasma' . 
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Fig. 1. — Background ID stratification for tliree models (Simulations 0, 1, 3). Dotted 
line: Simulation 0, the fully ionized simulation with /i^ = mi/2. Solid line: Simula- 
tion 1, the fully ionized simulation with = ^i- Dot-dashed line: Simulation 3, the 
partia lly ionized model. The red line is from a standard model of the solar convection 
z one dstix 2004). Th e red crosses are the ID semi-empirical atmospheric model of 



Vernazza et al. 



( Il98ll ). Panel (a): Temperature. Panel (b): Gas density. Panel (c) 



Reduced mass normalized by the proton mass. Panel (d): Gas pressure. 
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Fig. 2.— Panel (a): Initial po/B"^ at y = Utube for the initial magnetic flux tube 
embedded in the convection zone for the three small scale simulations, Simulation 1 
(the fully ionized simulation), Simulation 2 (the partially ionized simulation without 
ionization effects in the equation of state), and Simulation 3 (the complete partially 
ionized model). Panel (b): Initial density perturbation pi/po, as in Equation (IMl) . in 
the tube for the same three simulations. 
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Fig. 3. — Evolution of three small scale (ephemeral region size) flux tubes. The white 
lines show the in-plane field (constant contours of A^) and the color contours show 
the log of density. The contours are at seven values regularly spaced between 
(and inclusive of) the extrema of Az and 0.1 of this extrema (the extrema contour 
is a single point located at the center of the tube). Top row: Simulation 1, the 
fully ionized simulation. Middle row: Simulation 2, the partially ionized simulation 
without ionization effects in the equation of state. Bottom row: Simulation 3, the 
complete partially ionized model. 
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Fig. 4. — Panel (a): Height of the center and edge of the flux tubes as a function of 
time. The center is defined as the location of the minimum in the vector potential 
Az- The edge is defined as the intersection of Az = 0.1 min {Az) with the x = line. 
Panel (b): Rise speed in the convection zone of the centers of the flux tubes for the 
same three simulations. The curves for Simulations 1 and 2 lie on top of each other 
in panel (b). 




Height (Mm) Height (Mm) Height (Mm) 

Fig. 5. — Line plots along x = for the stabilizing (— 1/3(5) and destabilizing 
{—Hp-^ln(B)) terms in the magnetic buoyancy instability inequality (l35ll . Left: Sim- 
ulation 1, the fully ionized simulation. Center: Simulation 2, the partially ionized 
simulation without ionization effects in the equation of state. Right: Simulation 3, 
the complete partially ionized model. 
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Fig. 6. — Same as Figure [3] but at later times of 1457 s and 1967 s, showing later 
stage emergence of the flux tube in Simulations 1, 2 and 3. 
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Fig. 7. — Magnetic Reynolds number (using the Pedersen resistivity) as a function of 
height at X = for the two partially ionized simulations. The plot is taken at time 
1359 s for Simulation 2, and at time 1724 s for Simulation 3. 
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Fig. 8. — In-plane flux above 1.2 Mm as a function of time, normalized to the in-plane 
flux in the initial flux tube, for Simulations 1, 2 and 3. 
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Fig. 9. — Temperature as a function of height along the hne a; = 0. The times 
in each simulation are chosen so that the transition regions in each simulation are 
approximately co-spatial in height. 
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Fig. 10. — Panel (a): Change in mass above 1.2 Mm, normalized to the total mass in 
the initial flux tube, as function of time. Panel (b): Egh^ar, from Equation fl38p for 
the same three simulations. 




and Simulation 3 at 1920 s (panel c). The times are shown so that in each figure the 
0.1 min (Az) contour has reached approximately the same height. 
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Fig. 12. — Panel (a): Initial po/B"^ aty = ytube for the large scale flux tube simulations. 
Solid line is Simulation 4, the /3 = 4 fully ionized model. Dot-dash line is Simulation 
5, the /3 = 4 partially ionized model. Double-dot dash line is Simulation 6, the fully 
ionized model with /3 ~ 40. Long-dashed line is Simulation 7, the partially ionized 
model with f3 ^ 40. Panel (b): Initial density perturbation in the tube for the same 
three simulations. 
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(a) : Sim4 t= 1 578 s 



(b) : Sim4 t= 2307 s 



(c) : Sim5t= 1093 s 



(d) : Sim5t= 1821 s 
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Fig. 13. — Emergence of large scale (active region size) flux tubes, showing in-plane 
field (constant contours of A^) and color contours of the log of density. The 
contours are at seven values regularly spaced between (and inclusive of) the extrema 
of Az and 0.1 of this extrema (the extrema contour is a single point located at the 
center of the tube). Top row: Simulation 4, the fully ionized simulation with /3 = 4. 
Bottom row: Simulation 5, the partially ionized simulation with /3 = 4. 
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(a) :Sim6 t= 4615 s 



(b) : Sim6 t= 5829 s 
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(c) : Sim7 t= 3594 s 



(d) : Sim7 t= 4687 s 
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Fig. 14. — Emergence of large scale (active region size) flux tubes, showing in-plane 
field (constant contours of Az) and color contours of the log of density. Top row: 
Simulation 6, the fully ionized simulation with /3 = 40. Bottom row: Simulation 7, 
the partially ionized simulation /3 = 40. 
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Fig. 15. — Panel (a): Height of the center and edges of the flux tube as function of 
time. The curves for the tube centers in Simulations 4 and 5 lie on top of each other. 
Panel (b): Axial rise speeds in the convection zone for the same four simulations. 
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Fig. 16. — The amount of in-plane flux above 1.2 Mm as a function of time, normalized 
to the integral of below ?/ = at t = s, for Simulations 4, 5, 6 and 7. 




Fig. 17. — Temperature as a function of height at x = for all four large scale 
simulations. The times in each simulation are chosen so that the transition regions 
in each simulation are approximately co-spatial in the y-direction. 
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Fig. 18. — Panel (a): Change in mass above 1.2 Mm, normalized to the total mass 
in the initial flux tube, as a function of time. Solid line: Simulation 4, the fully 
ionized simulation with /3 = 4. Dot-dashed line: Simulation 5, the partially ionized 
simulation with the same /3 as Simulation 4. Double-dot-dashed line: Simulation 6, 
the fully ionized model with /3 = 40. Long dashed line: Simulation 7, the partially 
ionized simulation with the same /3 as Simulation 6. Panel (b): Eshear, from Equation 
( 138|) for the same four simulations. 



